*Note: this analysis was run using Stata 13.1 MP
use "survey data 4 replication.dta"

*Table A1 - Summary statistics
summ net_comm_support_01 church_attend_01 Catholicint_01  ///
red_priest_shr red_bishop5056 inst_vulner_ind maj_imper capital_invest_ind prod_heavyind prod_lightind prod_produce ///
red_priest_exp617 red_priest_exp18plus comm_support_01 solidarity officialsolidaritymembershiprate ///
primvoc_edu sec_edu hi_edu village_imp smtown_imp medtown_imp lgtown_imp city_imp occup_imp maleimp malemis age_corr month monthsq if net_comm_support_01~=.&church_attend_01~=., sep(0)

*Figure 1 --> see map files

*Figure 2 & Table A2 - Drivers of infiltration (wojew-level data)

*Figure 3 & Table A3 - Infiltration effects
xtmixed church_attend_01 red_priest_shr red_bishop5056 inst_vulner_ind i.maj_imper capital_invest_ind prod_heavyind prod_lightind prod_produce primvoc_edu sec_edu hi_edu village_imp smtown_imp medtown_imp lgtown_imp city_imp i.occup_imp  maleimp malemis age_corr  i.month || wojew_cbos:  || wojew_month: , pweight(wt_wojew)
outreg2 using taba3.doc, replace label dec(3) word symbol(**,*,#) nocons drop(capital_invest_ind prod_lightind prod_produce church_attend_01mis primvoc_edu sec_edu hi_edu smtown_imp medtown_imp lgtown_imp city_imp maleimp malemis age_corr  month monthsq)
margins, at(red_priest_shr=(.04)) saving(fileattend1) 
margins, at(red_priest_shr=(.19)) saving(fileattend2) 
margins, at(red_bishop5056=0) saving(fileattend3) 
margins, at(red_bishop5056=1) saving(fileattend4) 
margins, at(inst_vulner_ind =-.65)  saving(fileattend5)
margins, at(inst_vulner_ind =1.67)  saving(fileattend6)
margins, at(maj_imper = 1)  saving(fileattend7)
margins, at(maj_imper = 2)  saving(fileattend8)
margins, at(maj_imper = 3)  saving(fileattend9)

*Create space for graph
clear all
use fileattend1.dta
recode * (*=.)
save fileattend0.dta

*Make graph for labels
clear all
use fileattend1.dta
recode * (*=.)
replace _margin = 0 in 1
save fileattend00.dta

use "survey data 4 replication.dta"
combomarginsplot fileattend00 fileattend00 fileattend00 fileattend00 fileattend00 fileattend00 fileattend00 fileattend00 fileattend00 fileattend00 fileattend00 fileattend00, ///
	title ("Factor", color(gs15) size(medium)) ///
    labels("Russian Empire" "Prussian Empire" "Habsburg Empire" " " "High Organizational Vulnerability" "Low Organizational Vulnerability" " " "Patriot bishop" "No Patriot bishop" " " "High Patriot Priest share" "Low Patriot Priest share" ) ///
	recast(scatter) horiz ///
	ytitle(" ") yscale(range(0.5 9.5))  xtitle("I-I-I-I-I-I-I" "I-I-I-I-I-I-I", color(gs15) size(medium)) scheme(s2mono) 
	graph save labels.gph

combomarginsplot fileattend9 fileattend8 fileattend7 fileattend0 fileattend6 fileattend5 fileattend0 fileattend4 fileattend3 fileattend0 fileattend2 fileattend1 , ///
	title ("Church attendance", size(medium)) ///
    labels(" " " " " " " " " " " " " " " " " " " " " " " ") ///	
    recast(scatter) horiz ///
	ytitle(" ") yscale(range(0.5 9.5))  xtitle("Predicted church" "attendance", size(medium)) scheme(s2mono) 
	graph save attend.gph

	
xtmixed Catholicint_01 c.red_priest_shr##c.church_attend_01 red_bishop5056##c.church_attend_01 c.inst_vulner_ind##c.church_attend_01 i.maj_imper##c.church_attend_01 village_imp  capital_invest_ind prod_heavyind prod_lightind prod_produce primvoc_edu sec_edu hi_edu smtown_imp medtown_imp lgtown_imp city_imp i.occup_imp  maleimp malemis age_corr  i.month || wojew_cbos:  || wojew_month:church_attend_01 , pweight(wt_wojew)
outreg2 using taba3.doc, append label dec(3) word symbol(**,*,#) nocons drop(capital_invest_ind prod_lightind prod_produce church_attend_01mis primvoc_edu sec_edu hi_edu smtown_imp medtown_imp lgtown_imp city_imp maleimp malemis age_corr  month monthsq)
margins, dydx(church_attend_01) at(red_priest_shr=(.04)) saving(filetrust1) 
margins, dydx(church_attend_01) at(red_priest_shr=(.19)) saving(filetrust2) 
margins, dydx(church_attend_01) at(red_bishop5056=0) saving(filetrust3) 
margins, dydx(church_attend_01) at(red_bishop5056=1) saving(filetrust4) 
margins, dydx(church_attend_01) at(inst_vulner_ind =-.65)  saving(filetrust5)
margins, dydx(church_attend_01) at(inst_vulner_ind =1.67)  saving(filetrust6)
margins, dydx(church_attend_01) at(maj_imper = 1)  saving(filetrust7)
margins, dydx(church_attend_01) at(maj_imper = 2)  saving(filetrust8)
margins, dydx(church_attend_01) at(maj_imper = 3)  saving(filetrust9)

combomarginsplot filetrust9 filetrust8 filetrust7 fileattend0 filetrust6 filetrust5 fileattend0 filetrust4 filetrust3 fileattend0 filetrust2 filetrust1 , ///
	title ("Trust in Catholic Church", size(medium)) ///
    labels(" " " " " " " " " " " " " " " " " " " " " " " ") ///
	recast(scatter) horiz ///
	ytitle(" ") yscale(range(0.5 9.5))  xtitle("Church attendance effect" "on church trust", size(medium)) scheme(s2mono) 
	graph save trust.gph


xtmixed net_comm_support_01 c.red_priest_shr##c.church_attend_01 red_bishop5056##c.church_attend_01 i.maj_imper##c.church_attend_01 c.inst_vulner_ind##c.church_attend_01 capital_invest_ind prod_heavyind prod_lightind prod_produce primvoc_edu sec_edu hi_edu village_imp  smtown_imp medtown_imp lgtown_imp city_imp i.occup_imp  maleimp malemis age_corr  i.month  [weight=wt_comb] || wojew_cbos: || wojew_month:church_attend_01 , pweight(wt_wojew)
outreg2 using taba3.doc, append label dec(3) word symbol(**,*,#) nocons drop(capital_invest_ind prod_lightind prod_produce church_attend_01mis primvoc_edu sec_edu hi_edu smtown_imp medtown_imp lgtown_imp city_imp maleimp malemis age_corr  month monthsq)
margins, dydx(church_attend_01) at(red_priest_shr=(.04)) saving(filenetcomm1) 
margins, dydx(church_attend_01) at(red_priest_shr=(.19)) saving(filenetcomm2) 
margins, dydx(church_attend_01) at(red_bishop5056=0) saving(filenetcomm3) 
margins, dydx(church_attend_01) at(red_bishop5056=1) saving(filenetcomm4) 
margins, dydx(church_attend_01) at(inst_vulner_ind =-.65)  saving(filenetcomm5)
margins, dydx(church_attend_01) at(inst_vulner_ind =1.67)  saving(filenetcomm6)
margins, dydx(church_attend_01) at(maj_imper = 1)  saving(filenetcomm7)
margins, dydx(church_attend_01) at(maj_imper = 2)  saving(filenetcomm8)
margins, dydx(church_attend_01) at(maj_imper = 3)  saving(filenetcomm9)

combomarginsplot filenetcomm9 filenetcomm8 filenetcomm7 fileattend0 filenetcomm6 filenetcomm5 fileattend0 filenetcomm4 filenetcomm3 fileattend0 filenetcomm2 filenetcomm1 , ///
	title ("Net Communist Support", size(medium)) ///
    labels(" " " " " " " " " " " " " " " " " " " " " " " ") ///
	xline(0, lpattern(shortdash)) recast(scatter) horiz ///
	ytitle(" ") yscale(range(0.5 9.5))  xtitle("Church attendance effect" "on net communist support", size(medium)) scheme(s2mono)
	graph save commsupp.gph
graph combine labels.gph attend.gph trust.gph commsupp.gph, col(4) ti(x) xsize(8.5) title("Figure 3: Infiltration Effects vs. Alternative Explanations",  size(medsmall))  scheme(s2mono) ///
note("Note: Predicted values (left panel) and predicted marginal effects of church attendance (in center and right panel) with 95% confidence intervals." ///
"For full regression results see Table A3") 
graph save combined_effects.gph


*Figure 4 & Model 4 of Table A3 

xtmixed net_comm_support_01  c.red_priest_shr##c.church_attend_01##c.red_priest_exp617##c.red_priest_exp18plus red_bishop5056##c.church_attend_01##c.red_priest_exp617##c.red_priest_exp18plus i.maj_imper##c.church_attend_01 c.inst_vulner_ind##c.church_attend_01 capital_invest_ind prod_heavyind prod_lightind prod_produce primvoc_edu sec_edu hi_edu village_imp smtown_imp medtown_imp lgtown_imp city_imp i.occup_imp  maleimp malemis age_corr  i.month  [weight=wt_comb] || wojew_cbos:  || wojew_month: , pweight(wt_wojew)
outreg2 using taba3.doc, append label dec(3) word symbol(**,*,#) nocons drop(capital_invest_ind prod_lightind prod_produce church_attend_01mis primvoc_edu sec_edu hi_edu smtown_imp medtown_imp lgtown_imp city_imp maleimp malemis age_corr  month monthsq)
margins, dydx(church_attend_01) at(red_priest_shr=(.19) red_priest_exp617=0 red_priest_exp18plus=7) level(90) saving(fileexp6)
margins, dydx(church_attend_01) at(red_priest_shr=(.04) red_priest_exp617=0 red_priest_exp18plus=7) level(90) saving(fileexp5)
margins, dydx(church_attend_01) at(red_priest_shr=(.19) red_priest_exp617=7 red_priest_exp18plus=0) level(90) saving(fileexp4)
margins, dydx(church_attend_01) at(red_priest_shr=(.04) red_priest_exp617=7 red_priest_exp18plus=0) level(90) saving(fileexp3)
margins, dydx(church_attend_01) at(red_priest_shr=(.19) red_priest_exp617=0 red_priest_exp18plus=0) level(90) saving(fileexp2)
margins, dydx(church_attend_01) at(red_priest_shr=(.04) red_priest_exp617=0 red_priest_exp18plus=0) level(90) saving(fileexp1)


combomarginsplot fileexp2 fileexp1 fileexp4 fileexp3 fileexp6 fileexp5 , ///
	title ("Fig.4 Patriot priest exposure timing and church attendance effects", size(medium)) ///
    labels("No exposure - high infiltration" "No exposure - low infiltration" "Early exposure - high infiltration" "Early exposure - low infiltration" "Adult exposure - high infiltration" "Born Adult exposure - low infiltration", size(small)) ///
	xline(0, lpattern(shortdash)) recast(scatter) horiz ///
	ytitle(" ") yscale(range(0.5 6.5))  xtitle("Church attendance effect on net communist support", size(small)) scheme(s2mono) ///
	note("The figure shows the marginal effects and 90% confidence intervals of church" ///
		"attendance on the predicted change in the net communist support for different age cohorts ." )
graph save Graph "C:\Users\gpop\Google Drive\GreatDivide\submission 2020\JOPReviseResubmit\Figure 4 exposure diffs.gph"


*Table A4 - Disaggregated vulnerability results --> wojew-level data + do file


*Table A5 Response bias

meprobit solidarity c.red_priest_shr c.church_attend officialsolidaritymembershiprate capital_invest_ind prod_heavyind prod_lightind prod_produce church_attend primvoc_edu sec_edu hi_edu village_imp smtown_imp medtown_imp lgtown_imp city_imp i.occup_imp  maleimp malemis c.age_corr month monthsq  || wojew_cbos:  
outreg2 using tabrespbias.doc, replace label dec(3) word symbol(**,*,#) nocons drop(capital_invest_ind prod_lightind prod_produce church_attendmis primvoc_edu sec_edu hi_edu village_imp smtown_imp medtown_imp lgtown_imp city_imp maleimp malemis age_corr  month monthsq)

meprobit solidarity c.red_priest_shr##c.church_attend officialsolidaritymembershiprate capital_invest_ind prod_heavyind prod_lightind prod_produce church_attend primvoc_edu sec_edu hi_edu village_imp smtown_imp medtown_imp lgtown_imp city_imp i.occup_imp  maleimp malemis c.age_corr month monthsq  || wojew_cbos:  
outreg2 using tabrespbias.doc, append label dec(3) word symbol(**,*,#) nocons drop(capital_invest_ind prod_lightind prod_produce church_attendmis primvoc_edu sec_edu hi_edu village_imp smtown_imp medtown_imp lgtown_imp city_imp maleimp malemis age_corr  month monthsq)


*Table A6 - Borders --> wojew-level data + do file


*Table A7 Robustness tests

*Original specification (same as Table A3 model 3) - HLM with varying intercepts and slopes
xtmixed net_comm_support_01 c.red_priest_shr##c.church_attend_01 red_bishop5056##c.church_attend_01 i.maj_imper##c.church_attend_01 c.inst_vulner_ind##c.church_attend_01 capital_invest_ind prod_heavyind prod_lightind prod_produce primvoc_edu sec_edu hi_edu village_imp  smtown_imp medtown_imp lgtown_imp city_imp i.occup_imp  maleimp malemis age_corr  i.month  [weight=wt_comb] || wojew_cbos:  || wojew_month:church_attend_01 , pweight(wt_wojew)
outreg2 using tabrobust.doc, replace label dec(3) word symbol(**,*,#) nocons drop(capital_invest_ind prod_lightind prod_produce primvoc_edu sec_edu hi_edu village_imp smtown_imp medtown_imp lgtown_imp city_imp maleimp malemis age_corr  month)

*HLM with varying intercepts
xtmixed net_comm_support_01  c.red_priest_shr##c.church_attend_01 red_bishop5056##c.church_attend_01 i.maj_imper##c.church_attend_01 c.inst_vulner_ind##c.church_attend_01 capital_invest_ind prod_heavyind prod_lightind prod_produce primvoc_edu sec_edu hi_edu village_imp  smtown_imp medtown_imp lgtown_imp city_imp i.occup_imp  maleimp malemis age_corr  i.month  [weight=wt_comb] || wojew_cbos:  || wojew_month: , pweight(wt_wojew)
outreg2 using tabrobust.doc, append label dec(3) word symbol(**,*,#) nocons drop(capital_invest_ind prod_lightind prod_produce primvoc_edu sec_edu hi_edu village_imp smtown_imp medtown_imp lgtown_imp city_imp maleimp malemis age_corr  month)

*OLS with clustered standard errors
reg net_comm_support_01  c.red_priest_shr##c.church_attend_01 red_bishop5056##c.church_attend_01 i.maj_imper##c.church_attend_01 c.inst_vulner_ind##c.church_attend_01 capital_invest_ind prod_heavyind prod_lightind prod_produce primvoc_edu sec_edu hi_edu village_imp  smtown_imp medtown_imp lgtown_imp city_imp i.occup_imp  maleimp malemis age_corr  i.month  [weight=wt_comb_wojew], cluster(wojew_cbos)
outreg2 using tabrobust.doc, append label dec(3) word symbol(**,*,#) nocons drop(capital_invest_ind prod_lightind prod_produce primvoc_edu sec_edu hi_edu village_imp smtown_imp medtown_imp lgtown_imp city_imp maleimp malemis age_corr  month)

*Tobit with clustered standard errors
version 13
tobit net_comm_support_01  c.red_priest_shr##c.church_attend_01 red_bishop5056##c.church_attend_01 i.maj_imper##c.church_attend_01 c.inst_vulner_ind##c.church_attend_01 capital_invest_ind prod_heavyind prod_lightind prod_produce primvoc_edu sec_edu hi_edu village_imp  smtown_imp medtown_imp lgtown_imp city_imp i.occup_imp  maleimp malemis age_corr  i.month  [weight=wt_comb_wojew], ul ll cluster(wojew_cbos)
outreg2 using tabrobust.doc, append label dec(3) word symbol(**,*,#) nocons drop(capital_invest_ind prod_lightind prod_produce primvoc_edu sec_edu hi_edu village_imp smtown_imp medtown_imp lgtown_imp city_imp maleimp malemis age_corr  month)

*Using PZPR support instead of net communist support as DV 
xtmixed comm_support_01  c.red_priest_shr##c.church_attend_01 red_bishop5056##c.church_attend_01 i.maj_imper##c.church_attend_01 c.inst_vulner_ind##c.church_attend_01 capital_invest_ind prod_heavyind prod_lightind prod_produce primvoc_edu sec_edu hi_edu village_imp  smtown_imp medtown_imp lgtown_imp city_imp i.occup_imp  maleimp malemis age_corr  i.month  if net_comm_support_01~=. [weight=wt_comb] || wojew_cbos:  || wojew_month:church_attend_01 , pweight(wt_wojew)
outreg2 using tabrobust.doc, append label dec(3) word symbol(**,*,#) nocons drop(capital_invest_ind prod_lightind prod_produce primvoc_edu sec_edu hi_edu village_imp smtown_imp medtown_imp lgtown_imp city_imp maleimp malemis age_corr  month)
